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Summary. A rigid-plastic Cosserat model has been used to study dense, 
fully developed flow of granular materials through a vertical channel. Fric- 
tional models based on the classical continuum do not predict the occurrence 
of shear layers, at variance with experimental observations. This feature has 
been attributed to the absence of a material length scale in their constitutive 
equations. The present model incorporates such a material length scale by 
treating the granular material as a Cosserat continuum. Thus localised couple 
stresses exist and the stress tensor is asymmetric. The velocity profiles pre- 
dicted by the model are in close agreement with available experimental data. 
The predicted dependence of the shear layer thickness on the width of the 
channel is in reasonable agreement with data. In the limit of small e (ratio 
of the particle diameter to the half- width of the channel), the model predicts 
that the shear layer thickness scaled by the particle diameter grows as e — 1 / 3 . 



1. Introduction 

In many terrestrial flows of granular materials, gravity consolidates the medium 
to a state where sustained frictional contact between the particles is the dominant 
mode of momentum transfer. In this regime of high solids fraction and low de- 
formation rate, models based on concepts in metal plasticity and soil mechanics 
have been traditionally used to describe the flow Q . While many gross features of 
granular flows can be predicted using these models, one aspect they fail to capture 
is the thickness of shear layers; often when granular materials are sheared, large 
portions of the material do not suffer sustained deformation. In the experiments of 
Roscoe Nedderman and Laohakul ||, Gudehus and Tejchman Q, the velocity 
gradients are confined to layers approximately 5-40 particle diameters in thickness. 
Moreover, the thickness of the shear layers is influenced by the nature of the bound- 
aries; when the flowing medium is confined by smooth walls, it is found that the 
thickness of the shear layers is less than that in the case of rough walls ||, ||] . 

Conventional models of plasticity do not predict shear layers || f?J . The failure 
of the frictional models to predict the thickness of the shear layers accurately has 
been attributed to the absence of a material length scale in their constitutive equa- 
tions ||cited in D- To overcome this deficiency of the classical models, the particle 
size must be incorporated in the constitutive equations. In the absence of a com- 
prehensive micro-mechanical model to describe friction, a continuum theory that 
includes a material length scale in the constitutive equations can be constructed 
by modelling the granular material as a Cosserat continuum ||. We shall argue 
later that the frictional nature of particle interactions provides sufficient grounds 
for using this approach. We note here that models based on kinetic theory [see, for 
example, |l0|, involve the particle diameter in the constitutive relations. However, 
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these models are expected to hold only for rapid flows, where particle interactions 
may be approximated by instantaneous collisions. 

In this paper, we explore the use of a Cosserat plasticity model to describe 
steady, fully developed, plane flow of a granular material in a vertical channel 
under the action of gravity. The predictions of the model will be compared with 
data reported in the literature. While Cosserat plasticity models have been applied 
to problems in granular flow in the past ||, [| 0, O [l2| , these studies address 
unsteady flows and are posed in terms of strain increments; with this formulation, 
Tejchman and Gudehus |5j found it difficult to integrate the equations numerically 
for large times. To the best of our knowledge, the present work represents the first 
attempt to examine steady flow in this context. We indicate below how the model 
is developed, and then apply it to channel flow. 



2. The Cosserat model 

The field variables of the classical continuum are the density p, the linear velocity 
v, and the stress tensor a. A Cosserat continuum (jD|p. 223]; Q) involves two 
additional field variables, namely, the angular velocity uj, and the couple stress 
tensor M_. Considering a Cartesian coordinate system (figure 1), M xz represents 
the couple per unit area exerted about the z-axis on a plane x — constant, by 
the material to the left of this plane. A positive value of M xz is taken to impose 
an anti-clockwise rotation on this plane (figure 1). For a Cosserat continuum, 
the mass and linear momentum balances must be supplemented by the angular 
momentum balance, which relates w, A£, and a. For steady, fully developed flow, 
spatial gradients of M cause a to be asymmetric. This is in contrast to the classical 
continuum, which assumes implicitly that there are no couple stresses, body couples, 
and intrinsic angular momentum; hence the angular momentum balance can be 
satisfied identically by requiring a to be symmetric. 

There is enough analytical evidence to motivate the use of a Cosserat model in 
the present problem. Dahler [|l5| used a statistical mechanical approach to develop 
expressions for the stresses in a fluid composed of diatomic molecules. For molecules 
interacting via central forces, which are directed along the lines joining the centers 
of mass of the molecules, a is found to be symmetric. However, his model suggests 
that non-central forces may cause a to be asymmetric. Campbell [[l6) simulated 
the shearing of circular discs between parallel plates, assuming that the collisions 
between discs, and between a disc and the wall, were instantaneous. In the latter 
case, wall roughness was incorporated by imposing (after collision) either (i) a zero 
relative velocity between the surface of the disc and the wall, or (ii) a zero relative 
velocity between the center of the disc and the wall. In both cases a was asymmetric 
near the wall, and there were non-zero couple stresses. 

Jenkins et al. jl7j constructed a micro-mechanical model for an assembly of iden- 
tical spheres. They found that asymmetric stresses resulted when the distribution 
of contact normals was anisotropic; however, they secured the symmetry of the 
stress tensor by suitably enforcing the rotation of particles. 

Dry friction, the dominant mode of momentum transfer in high-density flows, 
introduces non-central forces in an inherently complex fashion. Hence, we expect 
that a micro-mechanical model for dry friction would result in a continuum with 



A FRICTIONAL COSSERAT MODEL FOR GRANULAR MATERIALS 



3 



asymmetric stresses; such materials can be modelled as Cosserat continua. A satis- 
factory micro-mechanical model is not yet available, and it is hoped that this issue 
will be addressed by future investigators. 

2.1. Equations of motion. It is instructive to write the equations for the case of 
steady plane flow, and later simplify them for the case of fully developed flow. For 
flow parallel to the xy plane (figure 1), the velocity field has the following form: 

v x = v x (x,y); v y = v y (x,y); v z = 0; u x = u y = 0; lo z = uj g {x,y), (1) 

where v x and lo x are the x components of the linear and angular velocity, respec- 
tively. A positive value of lo z is associated with an anti-clockwise rotation about 
the z-axis. 

The balances for mass and linear momentum are 

d_, . d_ 
dx dy ' 



ilK) + i:K) = ' ( 2 ) 



da xx d<7 yx ( d d . 

' p P v I v x — +v y — \ v x = 0, (3) 



dx dy \ dx dy 

d(7 X y d(Tyy 



dx • dy • P ^\ Vx d- X + ^d-y) V y = -P^ (4) 

where v is the solids fraction or the volume fraction of solids, cry's are the com- 
ponents of the Cauchy stress tensor, defined in the compressive sense, p p is the 
intrinsic density of the particles, assumed constant, and g is the acceleration due 
to gravity. 

Following Jaunzemis ]l3|,p. 233], the z component of the angular momentum 
balance is 

dM xz dM yz ( d d \ 

-q^T + —Qy- ^ Pv<* + - + Pv v + Vygjj ) Vz = 0, (5) 

where Mi Z s are the couple stresses, r\ z is the z component of the intrinsic angular 
momentum (per unit volume), and Q z is the z component of the body couple acting 
on the material. 

To close the above set of equations, constitutive relations for the cry and Mj Z 
are required. 

2.2. Constitutive equations. Miihlhaus and Vardoulakis || and Tejchman and 
Wu H have developed Cosserat plasticity models for studying the development of 
shear bands in granular flow. In their models, the yield condition and the flow 
rule were modified to account for the influence of the couple stress and to provide 
a relation for the angular velocity. We have adapted their model to the present 
problem. The constitutive equations comprise of a yield condition and a flow rule, 
which are elaborated below. 

2.2.1. Yield condition. Following Besdo (incited in [H, de Borst [|o), and Tejch- 
man and Wu || , we use a yield condition of the form 

F = t - Y = 0, (6) 

where 

1 \ 1/2 

WijVji + ( M Y ^3 M ijj J ( 7 ) 
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°ij = °w ~ (1/3) Cfcfe is the deviatoric stress, <5jj is the Kronecker delta, ai, 
a2, and L are material constants, and d p is the particle diameter. Here Ld p is a 
characteristic material length scale; the value of L will be chosen later, de Borst 



20 1 assumes that the yield limit Y depends on the mean stress a = <7fc/-/3, and a 
hardening parameter h. Here we identify h with the solids fraction v. 

Schoficld and Wroth |H],p. 112] and Jackson § discuss the use of a yield con- 
dition of the form F(q_, v) — in the classical frictional models. The yield con- 
dition (||) with Y = Y(a, v) represents an attempt to include couple stresses 
within this framework. Only two of the three parameters a\, a 2 and L in (]?]) 
are independent, because the third parameter may be absorbed in the definition of 
Y(a,v) (see (0)). Following Muhlhaus and Vardoulakis 0] we set a\ + a 2 — 1/2, 
without loss of generality. 

Tejchman and Wu || use A = aijax = 1/3 and L = 1. Here we retain their 
choice of A and and treat L as an adjustable parameter, whose value is chosen as 
described later, de Borst [^0) found that changes in the values of A and L affected 
the post-peak behaviour of a sample which was sheared between parallel plates. 
Unfortunately, neither experiments nor satisfactory micro-mechanical models are 
available to guide the choice of A. 

Following Prakash and Rao |2^] , we assume the following form for the yield limit 

Y 

Y = u c {y) sin <j> (na - (n - l)^"/^ 1 ^ ; a = — ^— . (8) 
V / a c (y) 

Here a c (y) is the mean stress at a critical state, is a material constant called the 
angle of internal friction, and n is a material constant. The significance of a critical 
state will be explained shortly. The dependence of a c on the solids fraction v is 
taken to be Q 

!0 v < u min , 

Av-^nY v . <v<v ( 9 ) 
T7i _ ,.\q ^mm v "max- 
V^max V ) 

Here A, v m \ n , f max , p and q are material constants. Note that a c [v) has been chosen 
to be zero below t'min, the solids fraction at loose random packing, and to diverge 
as v approaches ^ max , the solids fraction at dense random packing. 

2.2.2. Flow rule. Tejchman and Wu |(| have used incremental elasto-plastic consti- 
tutive equations, which they attribute to Muhlhaus [|lTJ . Elastic effects are ignored 
in the present work to simplify the analysis. Because we are interested in sustained 
flow, the plastic strain increments used by Tejchman and Wu Q are replaced by 
suitable velocity gradients |llj . In Cartesian tensor notation, the flow rule is given 

by 

dvi ■ dG dcui ■ dG 

or-' x n^r "■' or- x oJT/ (10) 

where G(a, M, v) is the plastic potential, e^-fc is the alternating tensor, and A is a 
scalar factor. We note here that is the sum of the rate of deformation tensor Dij 
and an objective antisymmetric tensor representing the difference between the spin 
tensor and the particle spin EijktUk- Eij and Hij are conjugate to the stress <jji and 
the couple stress Mji, respectively, in the sense that the rate of working per unit 
volume by the contact forces and couples is given by —(ajiEij + MjiHij) H]. In 
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a classical continuum, M vanishes and a is symmetric; hence the above expression 
reduces to —a^Dy, where Dij = (1/2) (dvi/dxj + dvj/dxi) denotes a component 
of the rate of deformation tensor. 

Lacking detailed information on the plastic potential G, we adopt the most 
commonly used closure, namely the associated flow rule [^l],p. 43]: 

G = F = t-Y. (11) 

This form of the flow rule accounts for density changes accompanying deformation. 



3. Application to channel flow 

For the case of steady, fully developed, plane flow, the velocity field is given by 

v y =v y (x); u> = u z (x), (12) 

and the other velocity components vanish. Hence the mass balance (||) is identically 
satisfied and the balances of linear and angular momentum reduce to 

-dx- = °> -dx- = - p ^ 9 > (13) 
dm „ ,\ 

-j^" + Vxy - <?yx = 0, (14) 

where m = M xz . It is assumed that the yield condition is satisfied at every point 
in the channel, so that the factor A in ([!(]) is always non-zero. In writing (Q), it is 
assumed that there are no body couples. 

3.1. The stress field. For fully developed flow, it will now be shown that all the 
normal stresses are equal. Equation ( |l0| ) implies that 

„ dv x A , , , , „ A dY , . 

E - = ^=°= 67^-^-^-3^ (15) 
Writing the corresponding equations for E yy and E zz and summing, we get 

dY 

or using (||) 

a = <j c (v)\ Y = <j c sm.(f). (17) 

Thus the material is at a critical state or a state of isochoric deformation, because 
En = V • v = 0. Comparison of (ra) and ( |l7|) shows that the value of n is not 
relevant. It also follows from (|l5|), ( |l7| ) and ^1^) that 

&xx = o yy = a zz = a c (v) = constant. (18) 

Hence the solids fraction does not vary across the width of the channel. 
Using ( p~7| ) and (|l8|), the yield condition (||) reduces to 

/ m 2 \ 1/2 

\ai{a 2 xy + <r 2 yx ) + 2a 2 a xy cjy X + — ^ 2 1 - a c {v) sin = 0. (19) 



6 



L. S. MOHAN, P. R. NOTT, AND K. K. RAO 



3.2. Velocity field. The non-trivial equations of the flow rule ( |l0| ) are 

E xy = oj = - (ai<j yx + aaVxy) , (20) 
dv y A 

E yx = -j uj = - (aicr X y + a 2 cryx) , (21) 

ax t 

H zx — -j— = — , , 2 ■ (22) 



On eliminating A we get 

dVy (A + l){a xy + <Ty x )u 

dx &yx A(7 ' X y 

du> 2(A+l)muj 



(23) 
(24) 



dx (Ld p ) 2 (a yx + Aa xy )' 

3.3. Boundary conditions. Considering symmetric solutions, we have 

a xy (0) = 0; w(0) = 0. (25) 

The angular velocity lo must vanish at the centerline of the channel, because a 
non-zero value implies a preferred direction of rotation. 
Equations (20) and ( p5| ) imply 

<j yx (0) = 0, (26) 

provided m(0) is bounded, i.e., v < ^ max (see (0), ([)[), and (|l9|)). Because a xy and 
<j yx both vanish at the centerline, the yield condition ( |l9| ) implies that the couple 
stress at the centerline is 

m(x = 0) = ±Ld p a c sin (j). (27) 

While both roots in (^) are mathematical solutions, only the negative root yields 
a physically reasonable solution. The justification for choosing this root, and the 
reason for discarding the other are discussed in Appendix A. 

At the right hand wall x — W we use the usual friction boundary condition ( [Q ; 

Hp- 40]) 

_£^=tan<5 at x = W, (28) 

where 5 is a constant called the angle of wall friction. Using (|l3| ) and (|l8|), J28|) 
reduces to 

^^tanJ at x = W, (29) 

which determines the value of v for specified values of W and 8. 
Following Tejchman and Gudehus 0], we assume that 

v y = —Kdpid at x — W 7 (30) 

where K is a dimensionless constant which reflects the roughness of the wall. To 
get a feel for this condition, consider a single spherical particle sliding or rolling 
down a vertical wall. Let v' y and cu' represent the linear velocity of the center of the 
particle and its angular velocity about an axis through its center, respectively. If 
the particle slides without rolling, lo' = 0, but v' y is arbitrary. Conversely, if it rolls 
without slipping \v' y \ = (d p /2) \lo'\. For the boundary condition (30), these limits 
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correspond to K — > oo and K — > 1/2, respectively. Reverting to the continuum, we 
expect that K will decrease as the wall roughness increases. 

One more boundary condition is needed to permit determination of all the inte- 
gration constants. Here we set 

u(x = W) = lu w , (31) 

where u w is a constant whose value is determined by adjusting either the mass 
flow rate or the centerline velocity to match the measured value. In experiments, 
the mass flow rate may be varied within limits by varying the width 2W S of the 
exit slot at the bottom of the channel. Because we are considering fully developed 
flow, W s does not occur explicitly in the governing equations, but its influence is 
incorporated by changing uj w . 

4. Solution procedure 
Introducing the dimensionless variables 

x d. 



£= 777; 



p . „, — v 



Tir\ 1/2 

W\ an m 



uj = u> \ — ; 0»i = 777; m 



9 J ' 11 P P gW' "" Pp gWd 

the balance equations (|l^) and ([l4]) may be rewritten as 

da 



' XX 



5 

doVi 



dm 

where v is the constant solids fraction across the channel, and 



0, (32) 

(33) 



e— + a xy - a yx = 0, (34) 



<=w (35) 



The dimensionless form of the yield condition (R3) is 



-2 



J 71 



ai(a xy + a yx ) + 2a 2 a xy a yx + —j = (a c sm(f>) . (36) 

Li 

Here a c {v) = a c /(p p gW). The flow rule (H|) and (|^) is given by 

du _ (1 + A) (a xy + a yx )U 



d£ a yx + Act xy 

du) 2(1 + A) 777 Hi 
di ~ eL 2 (a y 



The boundary conditions are: 
at the centerline (^ = 0) 



(37) 
(38) 



'xy 



0; <U = 0. (39) 
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at the wall (£ = 1) 

— j—r = tan 5; u — eKoJ. (40) 



Equations (26) and J27| ) may be written as 

W yx = 0, TFT — —NLv at £ = 0, (41) 

where N = sin0/tan<5, and ( f4(i| ) has been used to simplify the second of equa- 
tions @. 

4.1. Method of solution. For the special case A = —1, an analytical solution 
may be obtained as discussed in Appendix A. This shows that m = constant, 
~®yx = &xyi and uj and u — it(0) display a power law dependence on £. (The case 
A < — 1 is discussed in Appendix A.) 

We now discuss the case A > — 1. Inspection of (]32])-([4l"|) shows that the 
stress field is uncoupled from the velocity field. Hence we first integrate the equa- 
tions (|3^), ( |33| ) and (|3^). Equations (||^) and ( |33| ) may be solved along with the 
first of boundary conditions ( |39| ) to get 

~3xx — constant; a xy = — ^£. (42) 



The yield condition (|36j) may be solved for <j yx to get 

-2 \ \ V2 



AK : ( (A 2 - 1)K) 2 + 2(A + 1) ( N 2 v 2 - ^ ) j . (481 



In our calculations, only the root with the negative sign before the square root term 
in ([43]) was chosen. The justification for doing so is described in Appendix A. After 
substituting this expression for a yx in (|34|), we solve the equation along with the 
boundary conditions (|4l]) as an initial value problem by marching from £ = to 
£ = !• 

Equation (|34| ) is solved numerically using the LSODA routine p^j from ODEPACK 
in netlib. It should be noted that the above package estimates to at a small 
distance £1 from £ = as 

d.777 

« -NLv + — (0) £ a = —NLv, 
d£ 



because (41 ) and ( |34| ) imply that (0) = 0. This causes the term under the square 

root in (|]|) to be negative. To avoid this problem ([mJ) is integrated numerically 
from £ = £1 to £ = 1, with the initial condition given by 

with £1 = 10~ 5 . The use of a smaller value of £1 does not significantly affect the 
results. 



Here (0) is calculated by differentiating ( |34| ) with respect to £, and using (42) 

and (f43|). The resulting indeterminate expression is evaluated using the L'Hospital's 
rule to get 

e^J(0) = (A + l) v T -(A 2 - 1 + 2(A + 1)^^(0)) ^ (44) 
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This can be rearranged to get a quadratic equation for ^pr(O), and we choose the 
root that satisfies 

d 2 m , , 
e^r(O) < (A + 1)1/, 

as the other root is inconsistent with ([49|). 

Once the stresses are obtained, the velocities are calculated by integrating the 



flow rule (|37|)) and (g§) from £ = 1 to £ = using the initial conditions (flOj). The 



integration is started from £ = 1 because 



B(() = + ■ (45) 

£^ (.Cyx + Aa yx ) 



becomes unbounded as £ — > 0, and hence the right hand side of fl38| ) is indeterminate 
at £ = 0. 

Equation (B8h is therefore integrated from £ = 1 to get 



U(0 = ZJ W exp \- J B(0 d£'J , (46) 

where uJ w = uj w (W/g) 1 ^ 2 is the dimensionless angular velocity at the wall. It is 
shown in Appendix B that -B(£') d£' becomes unbounded as £ — > 0. Hence ZU 
satisfies the boundary condition ci7(0) = for all finite values of w w . 

It is also of interest to determine the behaviour of the solutions in the limit e = 
dp/W — > 0. The issue here is the scaling of the shear layer thickness as a function 
of the channel half-width, W. For small e, an asymptotic solution is constructed 
using a perturbation technique described in Appendix C. The predictions of this 
solution are discussed in the next section along with the numerical results. 

4.1.1. Parameter values. The intrinsic density of the particles (p p ) was taken from 
the studies of Nedderman and Laohakul ||, Natarajan et al. and Tuziin and 
Nedderman pg} ]. Glass beads were used in all the experiments. For want of data, 
the angle of internal friction <j> was taken to be equal to the reported angle of repose. 
The parameters v mm and v max were chosen to be 0.5 and 0.65, respectively. The 



parameters in (0) were estimated as follows. Jyotsna and Rao 29 used the data 
of Fickie et al. pG] to obtain an expression for the variation of the mean stress at 
a critical state (cr c ) as a function of the solids fraction v. This expression was used 
to generate the values of a c for v in the range 0.54-0.58, and the latter were used 
to estimate A/ (p p gW), p and q by the method of nonlinear least squares, using the 
Marquardt-Levenberg algorithm |^l|,p. 678]. The parameter values were found to 
be A/(p p gW) = 817, p = 2.5, and q = 2.2. 

In the experiments of Nedderman and Laohakul || and Natarajan et al. [27j| , a 
layer of particles was stuck to the walls of the channel. This will be referred to as 
a fully rough wall. When we compare model predictions with their data, the angle 



of wall friction is chosen as S = tan -1 (sin <j>) |32 . For comparing the predictions 
with stress measurements of Tiiziin and Nedderman p8| , the measured angle of 
wall friction, 5 = 10°, was used. 

The value of the parameter L, which occurs in the yield condition ( |36[ ) was esti- 
mated to be 10 by matching predicted velocity profiles with the data of Nedderman 
and Laohakul |3| (see figure 2) . This value was used in comparisons with all other 
data. The parameter K was set to 0.5. 
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5. Results 

In this section we compare the predictions of the theory with the data of Ned- 



derman and Laohakul 13], Natarajan et al. |27|, and Tiiziin and Nedderman |28 



5.1. Velocity profiles. With L = 10, there is a good match between predicted 
and measured linear velocity profiles of Nedderman and Laohakul || (figure 2). 
(Predictions of the theory with L = 2 and L = 20 are also shown in figure 2 
for comparison.) The solids fraction of 0.60 predicted by the model is in close 
agreement with the measured average value of 0.61. 

While there is no sharply defined plug layer in the model (and in the experi- 
ments) , there is a region of low shear rate near the center of the channel. In order 
to compare predictions with the data of Nedderman and Laohakul || , the apparent 
thickness of the "plug" layer, £ p is calculated from 

=0.95. (47) 

u(£ = 0) 

Hence the shear layer thickness, scaled by the particle diameter is A = (1 — 
£p)(W/d p ). The model predicts a central plug layer and a shear layer adjacent 
to the wall whose thickness is about 10.5 particle diameters. 

With L = 10, the predicted velocity profile also agrees well with the data 
of Natarajan et al. j2?j as shown in figure 3. This is an encouraging result be- 
cause the ratio of the channel width to the particle diameter differs by a factor of 
3.5 for the two sets of data. For the profile shown in figure 3, the solids fraction of 
0.59 lies in the range 0.55-0.67 estimated from the experiments. 

The open circles in figures 2 and 3 show the asymptotic velocity profiles for small 
e — the deviation from the numerical solution is greater in figure 3 because e is larger 
than that in figure 2. For e = 1/600, the asymptotic solution is indistinguishable 
from the numerical solution, as shown in figure 4. 

The angular velocity (uJ) profile, shown in figure 5, differs slightly from that of 
half the dimensionless vorticity (1/2) du/d£. As expected, the difference is more 
pronounced in the shear layer. (In a classical continuum, (l/2)du/d£ represents the 
local angular velocity of an infinitesimal spherical material volume.) The asymp- 
totic solution deviates significantly from the numerical solution for e = 1/30 (fig- 
ure 5), but the two solutions agree well for e = 1/600 (figure 6). 

5.2. Influence of channel width on the thickness of the shear layer. For a 

fixed value of the particle diameter d p , the thickness of the shear layer A increases 
with the half- width of the channel W (solid line in figure 7). This is roughly in 
accord with the data of Nedderman and Laohakul ||, which are represented by 
solid symbols in figure 7. For each value of the W/d p , there are three data points; 
these correspond to the estimates of A obtained by fitting three different functional 
forms to the measured velocity profile. 

For small values of e — d p /W, the perturbation solution (Appendix C) shows 
that 

fr 2x1/3 



Thus the dimensional thickness of the shear layer is proportional to (d p /W) 1 / 3 
when d p /W <C 1, and hence does not attain a constant value in this limit. It would 
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be interesting to conduct experiments with much larger values of W/d p than in the 
range shown in figure 7. This would permit a more stringent test of the model 
predictions. 

5.3. Influence of the parameter L on the thickness of the shear layer. 

As mentioned earlier, the length scale Ld p was chosen to fit the model predictions 
to the data of Nedderman and Laohakul It is important to know how the 
predictions vary with changes in this parameter. Figure 8 shows that the thickness 
of the shear layer is a weak function of L. In the limit of small e, the shear layer 
thickness varies as L 2 / 3 (Appendix C). 

5.4. Influence of the wall-roughness factor K on thickness of the shear 
layer. The variation of the shear layer thickness with the roughness parameter K 
is shown in figure 9. As mentioned earlier, K — > oo corresponds to a very smooth 
wall; it decreases as the wall roughness increases. Figure 9 shows that there is little 
variation with K of the shear layer thickness for small K, but significant variation 
in the range w 1-200. For K greater than 200, the velocity at the wall is greater 
than 95% of the centerline velocity. Hence, by our definition (fi"7j), the thickness 
of the shear layer is zero. As shown in Appendix C, the shear layer thickness is 
independent of K in the limit of small e. 



5.5. Stresses. 



5.5.1. Stress profiles. Figure 10 shows the profiles of the shear stresses ~d xy and ~5 yx 
for e = 1 /30 and 1 /600 . It is clear that the difference between a xy and W yx increases 
with £. Because a yx > ~d X yi the couple stress m also increases with £ (figure 11), 
in accord with ( 

The open symbols in figures 10 and 11 represent the asymptotic solution for small 
e. When e = 1/600, it is clear that the asymptotic solution is indistinguishable from 
the exact solution, and the difference a xy — a yx is also very small. 



34). 



5.5.2. Wall stresses. We now compare the predicted wall stresses with the data 
of Tiizun and Nedderman [g8[ (Table 1). The normal and shear stresses are over- 
predicted, but are of the same order of magnitude as the measured values. As noted 
by Mohan et al. (j , the dimensions of the channel used in the experiments are such 
that the front and the back faces may support a significant part of the weight of the 
material. Hence the shear stress measured at the side wall is expected to be less 
than the prediction, which assumes a channel of infinite depth. It is interesting to 
note that the estimate of Tiizun and Nedderman 1 28 for the average solids fraction 
is 0.63 and the model predicts a value of 0.625. 



6. Comparison with other models 

6.1. The classical frictional model. The classical frictional model predicts a 
flat velocity profile. This is consistent with the profile predicted by the Cosserat 
model in the limit d p — > for a fixed value of W. Further, the Cosserat continuum 
reduces to the classical continuum in this limit, because m — ► and a xy — > <r yx . 
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H Normal stress Shear stress 

Static Flowing Static Flowing 
0.91 2.7-3.6 4.0-5.1 0.30-0.39 0.29-0.35 
Experiments 1.07 2.7-2.9 3.8-4.6 0.30 0.41-0.58 
1.22 2.7-3.6 4.6-5.7 0.52-0.65 0.58-0.76 
Cosserat model 7.97 1.12 

Kinetic model 2.51 1.15 

Table 1. Comparison of predicted wall stresses with the data 
of Tiiziin and Nedderman Q for glass beads. Here H is the 
depth measured from the top of the channel. Units: H— m, 
stress- kN/m 2 . Parameter values: W = 0.155 m, d p = 2.29 mm, 
p p = 1180 kg/m 3 , (j> = 30°, 6 = 8°. 



6.2. The kinetic and frictional-kinetic models. The broken curves in fig- 
ure 7 show the results obtained by using the (high density) kinetic model and 
the frictional-kinetic model. For these models we have used the equations given 
in Mohan et al. [m] , except that the mean stress at critical state (a c in their paper) 
is evaluated using (||). 

The kinetic model is based on constitutive equations derived by using the kinetic 



theory of dense gases [see, for example, |10 . Two of the underlying assumptions of 



this theory, namely instantaneous binary collisions between particles and molecular 
chaos with respect to particle velocities, are expected to break down at high solids 
fractions. Therefore it is surprising, and perhaps fortuitous, that the predicted 
thickness of the shear layer is in fair agreement with the data (figure 7) even though 
the solids fraction is in the range of 0.64-0.65. 

Based on the results shown in figure 7, it is difficult to discriminate between the 
Cosserat and the kinetic models. It should be noted that both these models contain 
a material length scale in their constitutive equations. As noted by Tejchman and 
Wu [O], this may be a pre-requisite for a satisfactory description of shear layers. 

The frictional-kinetic model is constructed by assuming that the stress tensor 
is the sum of the kinetic stress tensor and the frictional stress tensor. This model 
grossly underestimates the thickness of the shear layer (see the dotted curve in 
figure 7) , probably because (i) frictional effects dominate kinetic effects in the shear 
layer, and (ii) the frictional constitutive equations do not contain a material length 
scale. 

6.3. The model of Tejchman and Gudehus [Q. The work of Tejchman and 
Gudehus §] appears to be the only other study which uses a Cosserat model for 
channel flow. They use an elasto-plastic model to examine the batch discharge of 
material from a cylindrical bin. The constitutive equations involve the Jaumann 
stress rate and the 'velocity strain' tensor. Since they do not present results for 
steady fully developed flow, a direct comparison of our predictions with theirs is not 
possible. We are currently attempting to use their model for the problem at hand, 
but some issues require consideration before results can be obtained. For example, 
Dienes [j33| have reported that the Jaumann stress rate furnishes an unrealistic 
oscillatory response in simple shear for a hypoelastic model. 
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7. Discussion 

Unlike the classical frictional model j7j, the present Cosserat model predicts 
velocity profiles which agree well with the data of Nedderman and Laohakul [|| 
and Natarajan et al. (^7). Further, the variation of the thickness of the shear layer 
with the width of the channel is also captured reasonably well by the model. The 
predicted wall stresses are of the same order of magnitude as the measured values, 
but there is considerable scope for improvement. In this context, it may be desirable 
to account for the finite spacing between the front and back walls. 

Our solution of the model for the limiting case of an infinitely wide channel (e — > 
0) with fully rough walls indicates that the shear layer thickness, scaled by the 
particle diameter, grows as e -1 / 3 , where e is the ratio of the particle diameter to 
the channel width. It would be interesting to compare this result with experiments 
conducted for a wide range of e. 

In the present work, and in most applications of the frictional Cosserat mod- 
els, ad hoc values are prescribed for the parameters a%, a,2, and L in the yield 
condition (0). Either suitably designed experiments, or a micro-mechanical treat- 
ment, would be valuable in providing estimates for these parameters. Similarly, it 
would be desirable to have a micro-mechanical basis for the kinematic boundary 
condition (p0|). An unsatisfactory feature of our model is that the solids fraction 
is constant across the channel. This is in variance with qualitative observations 
of Natarajan et al. |27j that the density in the shear layer is lower than that in the 
plug region. Perhaps the inclusion of elastic or kinetic effects in the model would 
correct this feature. In any case, accurate density measurements in channel flow 
are lacking, and more investigations in this direction are needed. 



Appendix A. Choice of signs in (f27J) and (|43|) 

In sections ^ and |I] it was mentioned that only the negative root of (27) was 
chosen in our calculations, and that only one of two possible choices was made in 
the sign for the square root in (^). In this appendix we discuss the justification 
for doing so. 

Rewriting the angular momentum balance ( |34| ) by substituting for a yx from (ff3|), 
we get 

e^ = {A + l)vt + D 1 '* = E+, (48) 
e |^ = (A + lK-£ 1/2 ^-, (49) 

where 

D = (A 2 - \)v 2 e + 2(A + 1) (n 2 u 2 - ^ , (50) 

and N = sin</>/ tan 6. Equations ( E^ ) and (fi"9| ) have to be integrated subject to the 
initial condition 

m(0) = ±NLv, (51) 



The qualitative behaviour of the solutions t o fl48 ), (|49|), and ( pip may be understood 
by examining the phase plane of (p8[) and fl4S|), such as that shown in figure 12. 



14 



L. S. MOHAN, P. R. NOTT, AND K. K. RAO 



For A < 1 - 2N 2 , D < at £ = 1. Hence @ and @ cannot be integrated till 
the wall (f = 1). 

For A > — 1 a real valued solution for ( fi"8| ) cannot be constructed using the 
initial condition to(0) = NLv, because E+ > in the region bounded by the 
curves D = 0, and D < outside this region (see figure 12). Suppose that the 
other initial condition m(0) = —NLv is used. For sufficiently small values of e, the 
trajectory touches the upper curve D — before reaching the channel wall (see, for 
example, the dot-dashed line in figure 12). Hence ( fl8| ) is discarded. Using a similar 
approach, it follows that a suitable solution can be constructed for ([l9j) only when 
the initial condition m(0) = —NLv is used. A typical trajectory is shown by the 
dotted line in figure 12. This choice of roots works only for small values of e; for 
the parameters used in figure 12, solutions could not be constructed for e = 0.33. 

For 1 — 2N 2 < A < — 1, and small enough e, it is possible to construct a solution 
for (fl8|), subject to the initial condition m(0) = NLv. 

For A = — 1, it follows from @ and @, and @ that m = constant = -NLv, 
W xy = W yx = —i/£, and the velocity profiles are given by 

57 = 57 l0 £ 2W /(« £ ), 

u = u(0) r, ; • 

V ' 2N/(eL) + 1 

The other root m(0) = NLv is discarded because lJ(0)/uj w — > oo as £ — > 0. 



Appendix B. Integration of 

In § 4.1 it was noted that 5J(0) = because lim£_>o B(€') d£' = oo, where 
B(£) is defined by (E^). This is shown below. 

Near £ = 0, the leading order behaviour of the stresses can be represented as 



m = m(0) + 0(£ 2 ), 



where 



da, 



u yx,0 - d< e 
The integral can now be written as 



5=o 



lim / d£' = lim 



" 2(A+l)m(0) 



,. 2(A+l)m(0). fa 
lim —7 — — u ' 



7 + / B(C)d£', 



In I j 



(52) 
(53) 



where a is a small positive number. When A > — 1, and 757(0) and a yx are evaluated 
as discussed in Appendix A, it follows from (|3^) and J44| ) that the factor multi- 
plying the logarithm in (^3|) is positive, and hence the expression in ( |53"| ) becomes 
unbounded. 
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Appendix C. Asymptotic solution for the Cosserat model 

Here we derive an asymptotic solution for the case of a fully rough wall (i.e., one 
for which tan<5 = sin0, or N = 1), in the limit e — > 0. 
Using (HI, (0), and (||), (§|) may be rewritten as 

- (A + IK + ((A 2 - 1) K) 2 + 2(A + 1) (V - |J) ) 1 12 = 0. (54) 

We now seek a solution for m of the form 

Til = mo + emi + e 2 t7i2 + • • • • (55) 



Substituting this in (|54|), expanding the term under the square root for small e, and 
collecting terms of 0(1) and 0(e) we get 

r 2 2 e 2 

m = -Lv{l-£ ) ; m x = - _ (56) 



Equation ( pq) shows that the solution for m is not uniformly valid as \emi/m^\ ~ 
0(1) when £ w 1 - (l/2)e 2/3 . 

To get an uniformly valid first approximation for m, we proceed as suggested by 
Van Dyke |34|p. 104]. Introducing new variables 



such that they are 0(1) in the inner region, we seek a solution of the form 

m' = m' + f(e)m[ + ... . (57) 

Substituting ( |57| ) in (|54"|), expanding the term under the square root for small e, 
and collecting terms of 0(1), we get 

dm'n _ „. m' 2 



£ = 2^ - (58) 



d£' vL 2 

Equation (|5^) is a Riccati equation |3^,p. 20], which can be converted to a second- 
order linear differential equation by using the transformation 

„ li(0 = ^W (59) 

to get 

d 2 m" 2£' „ 

5- = — 5-m . 60) 

d£' 2 £ 2 

Using the transformation £ = (2/L 2 ) 1 / 3 d60| ) reduces to the Airy equation 

d 2 TOn — ,, , ^ 

= Zml (61) 

and its general solution is given by |35],p. 100] 

mg = CiAi(|) + C 3 Bi(?). (62) 

Here Ai and Bi are the linearly independent Airy functions, and C\ and Ci are 
integration constants. Hence 

< = Of cAigffcBi,;) ? < C ' A ' (?) + aBi «» ■ < 63 > 
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To determine C± and C2, we follow the procedure discussed in Van Dyke |34|,p. 105]. 
The outer solution (|5^) is rewritten in terms of the inner variable and expanded for 
small e to get the leading order inner expansion of the outer solution 

Similarly, the inner solution should be rewritten in terms of the outer variables 
and expanded for small e. This is exactly equivalent to expanding ([33]) in the limit 
£ — > 00. The leading order outer expansion of the inner solution is 

-Lve 1 ^ (2£') 1/2 if C 2 = 0; 

Lvt 1 ^ (2£') 1/2 if C\ = 0, 

where the asymptotic expansions of the Airy function fl35|,p. 100] have been used. 
Thus the inner and outer expansions have the same functional behaviour in the 
"overlap" region provided C2 = 0; hence 

< = dAi (?). (64) 

Following the procedure described in Van Dyke |34],p. 94], the leading order com- 
posite (additive) solution is given by 

m = -Mi - ef' 2 + (§) + uLW - 0) 1/2 - (65) 

Similarly, the composite solution for the other variables is given by 

2e 2 \ 1/3 1 dAi 



= + 1 tt) j^-W' (66) 



2\ 1/3 



(67) 



u = 2lJ w [ — ) /(£), (68) 



where 



The parameter ui w = w(£ = 1) = cj(£ = 0) is determined as follows. Using ( p8| ) and 
the measured center line velocity u e (£ = 0), we get 

_ Ue(£ = 0) 

^ = 2X 1/3 ' (69) 

2e x 



'(Co) 

where £ = (2/L 2 ) 1 /3 £ -2/3, 

Using (p7|), d6S| ) and (|69|), the thickness £ p of the plug is given by 



where ? p = (l-£ P )(2/L 2 ) 1/3 e" 2 / 3 . 
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In the limit e — ► 0, £ — > oo and I(£ ) tends to a constant. Hence, for small e, 
£ p is approximately independent of e, and lim £ ^o£ P = 1-275. The dimensionless 
shear layer thickness is therefore given by 

fj 2xV3 



and the shear layer thickness expressed in terms of particle diameter is 

-1/3: 



W fL 2 ^ 1/3 



The above results are valid for the case of a fully rough wall (N — sin <fi/ tan S = 1). 
For smoother walls (N > 1), the outer solution is uniformly valid in the limit e — > 0. 
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Fig. 1. Elevation of the channel. The angular velocity lo z is pos- 
itive for an anti-clockwise rotation about the z-axis. The couple 
stress M xz exerted on the plane x = constant by the material to 
the left of this plane is positive when the couple is directed as 
shown. 
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Fig. 2. Scaled velocity profiles: (o) asymptotic solution, (■) data 
of Nedderman and Laohakul || for glass beads. The curves repre- 
sent numerical solutions with L = 2 ( ), 10 ( ) and 

20 ( ). Parameter values: A = 1/3, K = 0.5, W = 0.06 m, 

(e = 1/30), u(£ = 0) = 0.2, p p = 2940 kg/m 3 , and 
25° (5 = 22.91°). 



d p = 0.002 m 
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Fig. 3. Scaled velocity profiles: ( ) numerical solution, 

(o) asymptotic solution (■) data of Natarajan et al. for glass 
beads. Parameter values: W = 0.0255 m, d p = 0.003 m, (e = 0.12), 
dimensionless mass flow rate / wd£ = 0.15, cj) = 28° (S = 
25.15°), L = 10, the rest as in figure 2. 
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Fig. 4. Scaled velocity profiles: ( ) numerical solution, 

(o) asymptotic solution. Parameter values: W = 1.2 m, (e = 
1/600), L = 10, the rest as in figure 2. 
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f = x/W 

Fig. 5. Profiles of the angular velocity ( ) ZJ, 

( ) half the dimensionless vorticity, (l/2)du/d£, and 

the asymptotic solution for ZJ (o). Parameter values: L = 10, the 
rest as in figure 2. 
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Fig. 6. Profiles of the angular velocity: ( ) u>, and (o) the 

asymptotic solution for the angular velocity. Parameter values: 
W = 1.2 m, (e = 1/600), L = 10, the rest as in figure 2. 
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Fig. 7. Effect of the channel width 2W on the thickness of the 

shear layer (scaled by the particle diameter): ( ) numerical 

solution, (o) asymptotic solution, ( ) kinetic solution, and 

(- ) the frictional-kinetic solution. Here the solid symbols 

represent the estimates of Nedderman and Laohakul || , obtained 
by fitting the data to three different functional forms. Parameter 
values: L = 10, the rest as in figure 2. 
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Fig. 8. Influence of L on the predicted thickness of the shear layer: 

( ) numerical solution, (o) asymptotic solution. Parameter 

values as in figure 2. 




Fig. 9. Influence of K on the predicted thickness of the shear 
layer. Parameter values as in figure 2. 



28 



L. S. MOHAN, P. R. NOTT, AND K. K. RAO 




0.2 0.4 0.6 0.8 1 




Fig. 10. Profiles of the shear stresses: ( ) a xy , 

!jy X ( , e = 1/30; , e = 1/600). The symbols repre- 
sent the asymptotic solutions for a yx (o, e = 1/30; □, e = 1/600). 
Parameter values: L = 10, the rest as in figure 2. 
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Fig. 11. Profiles of the couple stress m ( , e = 1/30; 

e = 1/600). The symbols represent asymptotic solu- 
tions (o, e = 1/30; □, e = 1/600). Parameter values as in figure 2. 
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Fig. 12. Phase plane of the angular momentum balance equa- 
tions (Al) and (A2): , trajectory corresponding to (Al); 

, trajectory corresponding to (A2). In both cases, the ini- 
tial condition is m(0) = —NLv. Parameter values: N = 1, L = 10, 
(v = 0.603), A = 1/3, e = 1/30. 



